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Whole excitation spectrum is calculated within the Popov 
approximation of the Bogoliubov theory for a cylindrical sym- 
metric Bose-Einstein gas trapped radially by a harmonic po- 
tential. The full dispersion relation and its temperature de- 
pendence of the zero sound mode propagating along the ax- 
ial direction are evaluated in a self-consistent manner. The 
sound velocity is shown to depend not only on the peak den- 
sity, but also on the axial area density. Recent sound velocity 
experiment on Na atom gas is discussed in light of the present 
theory. 

PACS numbers: 03.75.Fi,05.30.Jp,67.40.Db,67.57.Jj 



Much attention has been focused on Bose-Einstein con- 
densation (BEC) since its realization in alkaline atomic 
gases. This macroscopic quantum phenomenon is unique 
because it occurs in ideal Bose gas, but real interest lies in 
the fact of how the particle interaction affects on it. Su- 
perfluidity which is absent in ideal Bose gas is understood 
to be sustained by the particle interaction. By flurrious 
experimental and theoretical studies, various fundamen- 
tal properties of Bose-Einstein condensation, including 
the condensation fraction relative to non-condensate, its 
temperature dependence and few lowest-lying collective 
modes have been already elucidated quite thoroughly [Q . 
However, because of the small limited size of trapped 
atomic gases these collective modes are all localized so 
far, thus the spectrum was discrete |^. 

Rapid development of experimental techniques to trap 
larger number of atoms leads to a chance to investigate 
propagating collective modes, such as zero, first and sec- 
ond sound. In particular, the zero sound mode with the 
continuous spectrum is intrigue because it gives the crit- 
ical superfluid velocity above which the system returns 
to be normal gas. In quantum liquid of superfluid ""^He 
where phonon, roton and maxon are known to exist mi- 
croscopic description for these modes has not been quite 
successful because "^He atoms are strongly interacting. In 
contrast, in the present dilute Bose gases there is a good 
opportunity to thoroughly understand these low-lying ex- 
citations in a microscopic level. In that context a recent 
zero sound measurement by Andrews et al. on a long 
cigar shaped Na atomic gas with length typically 500^m 
and width 10/xm provides an interesting testing ground 
and challenges theoretical investigation. 

The speed of this compressional zero sound mode is 
known to be given by c = \/nQg/m with the repulsive 
interaction constant g in the case of a homogeneous bulk 



system where the lowest collective mode is exhausted by 
the zero sound mode 0. Here no is the condensate den- 
sity. It is not obvious at all that this expression by re- 
garding riQ as the peak density is readily applicable to the 
above experiment as done by Andrews et al. B because 
the system trapped by a harmonic potential is finite and 
everywhere inhomogeneous. 

In fact, Zaremba S derives a different formula c = 
■y/ 710(0)5/2771 for an idealized cylindrical trap by solv- 
ing a time-dependent Gross-Pitaevskii equation within 
the Thomas- Fermi approximation (TFA) with 7io(0) be- 
ing the peak condensate density. The extra factor -\/2 in 
the denominator comes from a simple geometrical rea- 
son that the average condensate density is half the peak 
density in this parabolic shaped condensate. The same 
result is also obtained by Kavoulakis and Pethick [^. 
Zaremba's expression is far off the experimental data. 
As pointed out by himself a microscopic calculation is 
needed. 

Here we study the whole excitation spectrum, not only 
low-lying excitations, but also high energy ones for a 
cylindrical BEC system which is confined radially by a 
harmonic potential and infinite axially, approximating 
the above long cigar shaped trapped gas. We give spe- 
cial attention on zero sound velocity. The present set of 
the microscopic results, including condensate and non- 
condensate spatial profile, the relative fraction, their T- 
dependences and the whole excitation spectrum, within 
Popov approximation of the Bogoliubov theory enables 
us to draw a self-consistent picture of BEC systems. This 
theoretical framework has been quite successfully applied 
to the present dilute atomic BEC systems so far [y. 

We start out with the following Hamiltonian in which 
Bose particles interact with a two-body potential g5{ri — 
Y2) where g is a positive (repulsive) constant proportional 
to the s-wave scattering length a, namely g = AttU a/m: 



H= /"dr*^/i(r)^-|-| /"dr\^t^t^# 
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where the one-body Hamiltonian h{r) — „ 
/I, the chemical potential /i is introduced to fix the total 
particle number and V{r) is the confined potential. In 
order to describe the Bose condensation, we assume that 
the field operator ^ is decomposed into ^ = t/i -f where 
the ground state average is given by (^) — (j). A c-number 
(j) corresponds to the condensate wave function and 7/; 
is a q-number describing the non-condensate. Thus the 
condensate density is given by uq — |(/)p. The condition 
that the first order term in t/i vanish yields 



/i(r)0 + g|0|'<^ + 2gp(/) = O. 



(2) 



When the non-condensate density p is made to zero, it 
reduces to the Gross-Pitaevskii (GP) equation: h{r)(f) + 
g\4>\'^(t) = 0. The condition that the Hamihonian be di- 
agonalized gives rise to the foUowing set of eigenvalue 
equations for Uq and Vq with the eigenvalue Sq: 

{h{r) + 2gp + 2g\(j)\^}uq - gcb^Vq = SqUq (3) 

{h{r) + 2gp + 2g\(j)\^}vq - g^*^Uq = -EqVq. (4) 

The variational parameter p is determined self- 

COnsistently hy p = {ij^'^'lp) = J2q{u*qUq + VqVq)f{eq) + 

'^qVqVq whcrc /(s) is the Bose distribution function. 
These constitute a complete set of the equations within 
the Popov approximation. The total particle number 
density is given as rit — uq + p. 

We now consider a cylindrically symmetric system in 
the cylindrical coordinate: r = {r,9,z). The system 
is trapped by a harmonic potential V{r) — ^mojor^ 
radially and periodic along the z-axis whose period is 
L. The eigenfunctions of Uq and Vq are written as 



Uq{r)e 



no ^ikz 



and Vn 



Vq{r)e 



iw^rkz^ The set of the 



quantum numbers q is described by (k, Z, k) where k — 
0,1,2,---, / = 0,±1,±2,---, k = 0,±27r/i,±47r/L,---. 
Then, the functions Uq{r) and Vq{r) are expanded in 

terms of <^|''(r) = ^w\ Ji ("l*''"^) ^ "'?('') = 

Si Cg VJ (r) and Vq{r) — X^i '^9 V/ (^) where Ji{r) is 
the Bessel function of l-th order, aj denotes i-th zero 
of Ji, and i? is a sufficiently large radius. The eigenvalue 
problem reduces to diagonalizing a finite matrix. Note 
that since the angular momentum I and the wave num- 
ber k are both good quantum numbers, the eigenvalue 
equation is decomposed into each I and k. Each block- 
diagonal eigenvalue equation gives rise to the quantum 
number k along the radial direction. The iterative calcu- 
lations yield a converged self-consistent solution of Eqs. 
(2), (3) and (4) §. ^ 

The zero sound experiment by Andrews et al. ra is 
done under the fixed total number of Na atoms {'^ 
5 X 10®) and fixed axial trapping frequency {^ 18Hz) 
by varying the radial trapping frequency wq (typically 
~ 200Hz), giving rise to the zero sound speed as a func- 
tion of the peak condensate density no(0). To simulate 
this experiment we have performed self-consistent calcu- 
lations for various uj^ values and various area density n^ 
per unit length along the z-axis, which is uniform in our 
setting (for example, riz — 2 y. 10^//im). The scattering 
strength is chosen as a = 2.75nm Q. The area density 
Uz turns out to be one of the key parameters to fully 
characterize the experimental results, in particular, the 
zero sound velocity. The energy is scaled by wq- 
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FIG. 1. Spatial profiles of the condensate no(r) (full line), 

non-condensate p[r) (dotted line) and total nt{r) (dashed 

line) along the radial direction, n^ = 2 x 10* //^m, ujq = 200Hz 

and T = 0.6/xK. The inset shows the T-dependences of no(0), 

p(0) andnt(O). 

In Fig. |l| we show a typical example of the spatial pro- 
files of condensate no(r), the non-condensate p{r) and 
the total rit (r) along the radial direction at a finite tem- 
perature. The spatial extension of the gas is an order of 
10/im, which coincides with the typical radial width of 
a long cigar in the experiment p]. It is also seen that 
the condensate concentrates at the center r = 0, whose 
spatial variation is quite parabolic as predicted by TEA 
S while the non-condensate which is absent in TEA is 
pushed outward. The total density profile ritir) devi- 
ates far from a parabolic shape. The peak condensate 
density no(0) decreases progressively as T increases as 
shown in inset of Eig. Il|. These results are generally con- 
sistent with Hutchinson et al. for a spherical case who 
give detailed T dependence of these quantities within the 
same approximation. 
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FIG. 2. Excitation spectra (a) along the angular momen- 
tum I where the dotted line is the TFA prediction and (b) 
along the axial wave number k where the dotted line is the 
slope of the zero sound mode, riz — 2x W^/pm, ojq ~ 200Hz, 
and T = CB/iK. 

The whole excitation spectra characterized by the 
three quantum numbers q — {k, I, k) are shown in 
Eigs. 11(a) and (b). In Eig. ||(a) the energy levels are 
plotted as a function of the angular momentum I. The 
behavior e — uJovl expected by hydrodynamic theory (or 
TEA) due to Stringari ||] is only approximately obeyed 



when I is small in the present situation as shown by the 
dotted hne. 

The lowest three modes with Z = 0, 1, 2 for k = are of- 
ten discussed theoretically and observed experimentally 
on Rb case under a different geometrical confinement 
[^. The lowest dipole mode / — 1 coincides with loq, 
almost independent of the interaction strength because 
it corresponds to a rigid oscillation of the entire trapped 
gas. Our result satisfies this criterion as inspected from 
Fig. 0(a), ensuring the accuracy of our numerical cal- 
culations. Here we only remark that the relative posi- 
tion of the monopole mode e(0, 0, 0) with / = and the 
quadrapole mode e(0, 2,0) with / = 2 can be reversed 
when the area density along the z-axis changes where 
the energy eigenvalue is denoted by e(K, /, fc) in terms of 
the above quantum numbers. 

The zero sound velocity c is given by the slope of the 
lower edge of the spectrum along the k direction at fc = 
shown in Fig. ^(b). It is clear that c{k) is dispersive 
from a linear to quadratic curve, implying that the mode 
changes from the collective zero-sound to individual par- 
ticle excitations toward short wave numbers (Also see 
Fig. 0). This crossover point corresponds to the inverse 
of the characteristic coherent length ^ of the system. 
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FIG. 3. Zero sound velocity c against the peak conden- 
sate density no(0) at low T. The full curves are the har- 
monic potential case: tIz = 2 x 10 , 1 x 10 and 2 x 10 //im 
from top to bottom. The topmost one nearly corresponds 
to Zaremba's result. The dashed line is the anharmonic case 
(nz = 2 X 10'' //im). The dots are the experimental data H. 
The diamonds are the rigid wall potential case corresponding 
to infinite homogeneous system (dotted line). 

In Fig. ^ we plot the sound velocity c as a function 
of the peak density no(0) of the condensate. It is found 
that c depends on the area density Uz of the conden- 
sate along the axial direction, not only on no (0) . This is 
physically reasonable because along each curve in Fig. g 
the radial harmonic potential with ajp varies, which is 
implicit in this figure, keeping the area density n^ fixed. 
Therefore to fully characterize the system we need the 
additional parameter n^. In other words, c increases by 
increasing loq while fixing no(0). These curves should 
be compared with that in the infinite homogeneous case. 
The latter is calculated by the same scheme here except 
that a rigid wall boundary condition is employed instead 



of a harmonic potential. Our calculation indeed recov- 
ers the formula c = yriog/™ when the treated system 
is large enough where the all relevant quantities are spa- 
tially uniform except for the boundary M. In this case 
there is no additional parameter like n^. The velocity c 
increases as the area density n^ decreases as seen from 
Fig. 0. In the low area density limit c tends to approach 
Zaremba's result; c = ^JnQg/2ra. Therefore, within the 
present situation of the radial harmonic confinement our 
velocity never exceeds Zaremba's one. 

In Fig. ^ the T-dependence of c is displayed. As T 
decreases c(fc — 0) decreases, signaling the loss of the 
rigidity of this macroscopic quantum state toward the 
transition temperature. In the inset of Fig. Q the full 
dispersion relation c(fc) = ^ of the zero sound mode 
is displayed. The inverse of the crossover point from a 
linear to quadratic curve, that is, the departure point 
from a constant c{k) gives rise to the coherent length ^ 
(~ 0.8/im calculated for the present case). It is known 
that the coherent length is given by ^o 
homogeneous system. 
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FIG. 4. T dependence of the sound velocity. 

Tiz = 2 X 10''//^m and ujq = 200Hz. The inset shows the full 
dispersion relation c{k) of the zero sound mode as a function 
of fc at T= 1.2/iK. 

Let us now consider the experimental data for the zero 
sound by Andrews et al. [^ . For the reasonable area den- 
sity riz ^ O(10^)//xm our theoretical curves are far off the 
experimental data as seen from Fig. y. The experimental 
data are explained neither by our theoretical curves nor 
by the Zaremba curve based on TFA gH . We argue the 
possible origins of this discrepancy: (1) Uncertainty of 
the value of the repulsive interaction g or the s-wave scat- 
tering length a. We have assumed a — 2.75nm according 
to B . As inferred from the expression c^ oc a for infinite 
system, c is sensitive to the choice of the value a. If we 
choose a = 4.9nm as adopted before by MIT group pO[ , 
this discrepancy may be reduced. (2) Anharmonicity of 
the trapping potential. The possible anharmonicity such 
as kr^ in the trapping potential V{r) may be present. 
An example of our result for typically fc = 5 x 10~® J/m''^ 
is shown in Fig. ^ (dashed line), giving rise to impressive 
improvement of our fitting. Any anharmonicity of this 
type makes the condensate profile a flat-topped bulky 



shape, leading the system to more infinite homogeneous 
hke. This results in the more bulk like sound velocity. 
Note that the reduction factor \f2 in Zaremba's result 
comes from a simple geometrical reason characteristic to 
the parabolic shape condensate profile as mentioned be- 
fore. (3) As is seen from Figs. and the zero sound 
mode is dispersive. The true zero sound measurement 
is only realized when exciting long wave length limiting 
mode. Any deviation from this limit tends to increase 
the estimate for the velocity as seen from inset of Fig. H. 
This might be one reason why the low density data de- 
viates even from the infinite case, which gives the upper 
limit velocity. Here we point out that since the low den- 
sity case has smaller aspect ratio, the sound propagation 
might not be purely axial mode, possibly exciting the ra- 
dial mode simultaneously |lj]. (4) The actual system is 
confined axially even though it is very long. Thus along 
the axial direction the density is varying spatially, this 
making experiment and theory for the zero sound prob- 
lematic. We have performed the additional computation 
to see the shape and the area density in the actual exper- 
imental situations where the gas is trapped axially by a 
weak harmonic potential. The Gross-Pitaevskii equation 
mentioned before is solved under the axial trapping fre- 
quency 18Hz fixed for various radial trapping frequencies 
Wo as in the actual experiment [p[. From Fig. B which 
shows typical examples of the condensate spatial shape, 
it is seen that upon varying loq not only the peak density 
7io(0) but also the length and width of the condensate 
change, resulting in the change of the area density n^. 
Generally, n^ increases as wqi or no(0) increase. Thus 
in the actual experiment the velocity data are traversing 
the various curves in Fig. toward the upper ones as 
no(0) increases. We also note that according to Isl the 
sound velocity is position-dependent and increases near 
the end of the cigar. This is understandable because the 
"effective" density consisting of the actual density and 
the confining potential increases there, making c larger. 

In conclusion, we have calculated the whole excitation 
spectrum of a cylindrical symmetric Bose-Einstein gas 
trapped by an axial harmonic potential within the Popov 
approximation of the Bogoliubov framework. The zero 
sound compressional mode is analyzed in detail. The 
sound velocity and its temperature dependence are eval- 
uated. The recent experiment by Andrews et al. B on 
Na gas is discussed in light of the present theory, pointing 
out several possible origins of the disagreement between 
our theory and their result on the sound velocity. It is 
our hope that the whole set of physical measurable quan- 
tities determined here in a self-consistent manner further 
serves to experimentally check the internal consistency 
of the Bogoliubov framework, which has been remark- 
ably successful so far. We don't think the sound velocity 
an exception. 
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FIG. 5. Stereographic view of the condensate no(r, 2) at 
T — calculated by Gross-Pitaevskii equation at ujo = 200Hz 
(upper) and ujo = 600Hz (lower) under the axial potential 
fixed (18Hz). The total number is 5 x 10®. The area density 
riz at the center is seen to be 4.2xl0*/juni and 2.7xl0*/^ni 
respectively. 
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